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The properties of hydrogen-helium mixtures at high pressure are crucial to 
address important questions about the interior of Giant planets e.g. whether 
Jupiter has a rocky core and did it emerge via core accretion? Using path 
integral Monte Carlo simulations, we study the properties of these mixtures 
as a function of temperature, density and composition. The equation of 
state is calculated and compared to chemical models. We probe the accuracy 
of the ideal mixing approximation commonly used in such models. Finally, 
we discuss the structure of the liquid in terms of pair correlation functions. 

PACS numbers: 62.50.+p, 02.70. Lq, 64.30. +t 



1. INTRODUCTION 

Hydrogen and helium are the two most abundant elements in giant 
planets. While Jupiter and Saturn are well characterized on the surface, 
many basic questions about its interior have not been answered^. Jupiter's 
surface composition has been measured in situ by the Galileo probeP H 
74.2% by weight, He 23.1%, and 0.027% heavier elements, which is enhanced 
compared to the protosolar composition: 0.015% heavier elements along with 
H 73.6% and He 24.9%. The abundance of heavy element in the interior and 
their distribution are not well characterized.^ In particular, it is conversial 
whether Jupiter has a rocky core. The detection of a core in Jupiter may 
validate the standard model of giant-planet formation, nucleated capture 
of nebular hydrogenP An alternative scenario was proposed by Alan Boss 
who suggested that giant planets form directly from spiral instabilities in 
protostellar disks. Under this gravitational instability hypothesis, giant 
planets would not have a core, or at least a much smaller one. 

Since there is no direct way to detect a core in Jupiter, one must instead 
refer to indirect measurements and to models for the planet interior. Such 
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models are constrained by the available observation data, in particular, the 
properties at the planet surface and the gravitational moments measured 
through fly-by trajectories. All models rely on an equation of state (EOS) 
of hydrogen-helium mixtures. However, the uncertainties in the available 
EOS are large and have not allowed one, among many other questions, to 
determine whether Jupiter has a core. 

In this article, we present results from path integral Monte Carlo (PIMC) 
simulations^ that enable us to study quantum many-body systems at finite 
temperature from first principles. In this simulation, hydrogen-helium mix- 
tures are represented by an ensemble of electrons, protons and helium nuclei, 
each described by a path in imaginary time to incorporate quantum effects. 
The electrons are treated as fermions while exchange effects for the nuclei 
can be neglected for the considered thermodynamics conditions. 



2. PATH INTEGRAL MONTE CARLO 

The thermodynamic properties of a many-body quantum system at fi- 
nite temperature can be computed by averaging over the density matrix, 
p = e~P R ,/3 = l/k^T. Path integral formalisrrPis based on the identity, 
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where M is a positive integer. Insertion of complete sets of states between 
the M factors leads to the usual imaginary time path integral formulation, 
written here in real space, 

p(R, R'; 0) = J ... J 'dRi . . . tflW-i p(R, Ri; r) . . . p(R M -i, R'; r) (2) 

where r = (3/M is the time step. Each of the M steps in the path now has 
a high temperature density matrix p(Rfe,Rfe + i;r) associated with it. The 
integrals are evaluated by Monte Carlo methods. The density matrix for 
bosonic and fermionic systems can be obtained by projecting out states of 
corresponding symmetry. In PIMC, one sums up different permutations P, 

Pb/f (R,R;/3) = (±lfp(^VR';P) = ^£ (±lf J dR t e~ u ™. 

(3) 

For bosons, this is essentially an exact numerical algorithm that has found 
many applications.^ For fermions, the cancellation of positive and negative 
contributions leads to numerically unstable methods, which is known as the 
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fermion sign problem. Ceperley showed that this problem can be solved by 
introducing the restricted path approximationpE3 



MRo,R';/3)~^£(-ir / dEUe-oWl , (4) 
v 
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where one only samples path R(t) that stay within the positive region of 
a trial density matrix, pr(R(t), Ro; t) > 0. This procedure leads to an 
efficient algorithm for fermionic systems. All negative contributions to di- 
agonal matrix elements are eliminated.^ Contrary to the bosonic case, the 
algorithm is no longer exact since it now depends on the approximations 
for the trial density matrix. However, the method has worked very well in 
many applications! 12 * 1 ^ For pr, one can use the free particle or a variational 
density matrixP^In the results presented here, the high temperature density 
matrix was taken as a product of exact pair density matrices, 
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where po(R, R'; r) is the free particle density matrix and u(jCij, •; r) is the 
pair action for paths initially separated by r^j and finally at time r by r • • . An 
approximation is introduced by assuming that the different pair interactions 
can be averaged by independent Brownian random walks that are denoted 
by brackets (■■■)• This approach is efficient but not exact, and therefore 
puts a limit on the imaginary time step r in many-body simulations. The 
pair action, u, can be computed by different methods ! 15 * 1 ^ 17 ! The following 
simulation results were derived with N e = 32 electrons and the corresponding 
number of protons, N p , and helium nuclei Ane to obtain a neutral system 
(Ae = N p + 2A^Hc) with helium fraction x = 2N~n e /N e . Periodic boundary 
conditions are applied. As imaginary time discretization, we employ r = 
0.079. (Atomic units of Bohr radii and Hartrees will be used throughout 
this article.) 

The electrons are treated as fermions with fixed spin. We use variational 
nodes to restrict the paths and have therefore extended the approach in 
Ref. ^1] to mixtures. The nuclei are also treated as paths but their exchange 
effects are not relevant here. 
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Fig. 1. Density-temperature phase diagram of hot dense hydrogen. The blue 
dash-dotted lines separate the molecular, the atomic, the metallic, and the 
plasma regime. The green solid lines are isentropes for Jupiter and stars 
with 0.3, 1, and 15 solar masses. Single shock Hugoniot states as well as the 
inertial confinement fusion patrP^ are indicated by dashed lines. The thin 
solid line show p-T conditions of PIMC simulations. 



3. PHASE DIAGRAM OF HOT DENSE HYDROGEN 

Figure ^ shows the high temperature phase diagram of dense hydro- 
gen beginning with the fluid and reaching up to a highly ionized plasma 
state. The figure includes the isentropes for Jupiter and low mass stars^ 
and indicates the thermodynamic conditions, at which PIMC simulations 
have been applied! 2 ™ 21 ! 1 ^ We are now going to use these simulation results 
to characterize hot dense hydrogen from a path integral perspective. 

At low density and temperature, one finds a fluid of interacting hydrogen 
molecules. A PIMC snapshot for T = 5 000 K and r s = 4.0 is shown in Fig.|2j 
The proton paths are very localized due to their high mass. Their spread can 
be estimated from the de Broglie thermal wavelength, A 2 , = h 2 (3/2irm. The 
electron paths are more spread out but they are localized to some extent since 
two electrons of opposite spin establish the chemical bond in the hydrogen 
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Fig. 2. Snapshot from PIMC simulations of pure hydrogen in the molecular 
regime at T = 5 000K and r s = 4.0. The pink spheres denote the protons. 
The bonds (white lines) were added as a guide to the eye. The electron 
paths are shown in red and blue [light and dark gray] depending on their 
spin state. 

molecule. 

If the temperature is raised from 5000 to 10 6 K, hydrogen undergoes 
a smooth transition from a molecular fluid though an atomic regime and 
finally to a two-component plasma of interacting electrons and protons. 
Many-body simulations at even higher temperatures are not needed since 
analytical methods like the Debye-Hiickel theory work well. PIMC with ex- 
plicit treatment of the electrons can also be applied to temperatures below 
5 000K but groundstate methods are then more practical since excitation 
become less relevant and the computational cost scale with the length of 
paths like T _1 . 

A detailed analysis of the chemical species present in the low density 
regime (r s > 2.6) is given in Ref. 1201 With decreasing density, one finds 
that the degree of molecular dissociation increases since the atomic state 
has higher entropy. For the same reason, one observes that the degree of 
atomic ionization increases with decreasing density. All these low-density 
effects can be well characterized by analytical models based on approximate 
free energy expressions for atoms, molecules and ionized 
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Fig. 3. Deuterium in the dense molecular regime (T = 5 000K and r$ = 
1.86) Due to the density increase compared to Fig. |2] (see details there), 
the electron paths permute with a rising probability (shown as yellow [light 
gray] lines) but are still localized enough to form a bond between the two 
protons in the molecule. The electron density average over many electron 
configurations is indicated in gray color on the blue rectangles. 

In this regime, one also finds reasonably good agreement between for the 
EOS derived from PIMC and chemical modelsP^ 

If the density is increased at T = 5 000 K, one finds an intermediate 
regime of strongly interacting molecules (Fig. |3J). Some electron paths ex- 
change with neighboring molecules indicating the importance of fermionic 
effects. However, the electrons are still localized enough to provide a suffi- 
cient binding force for the protons. 

If the density is increased further from r s = 1.86 to 1.60, this binding 
force is lost due to further derealization of the electrons (Fig.HJ. Almost all 
electrons are now involved in long exchange cycles indicating a highly con- 
ducting, metallic state. No binding forces of the protons can be observed. 
One can therefore conclude that hydrogen metallizes in dissociated, or atomic 
form, a conclusion that is consistent with DFT simulation j 2 fr l 2 7 l 2 $ l 2 !fl if the 
density is increased further, the electrons form a rigid neutralizing back- 
ground and one recovers the limit of a one-component plasma of protons. If 
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Fig. 4. Deuterium in the metallic regime characterized by unpaired protons 
and a gas of degenerate electrons. The snapshot was taken from a PIMC 
simulation at T = 6 250K and r s = 1.60. The electron paths are delocalized 
and exchange frequently (see description of Fig. |3J) . 



the temperature is increased, the electron paths get shorter and shorter, the 
degree of electron degeneracy decreases gradually and one recovers the limit 
of a two-component plasma at high temperature. 

The nature of the molecular-metallic transition has not yet been deter- 
mined with certainty. A large number of analytical modelJ^ predict a first 
order transition, others do not. 23 PIMC simulation with free particle nodes 
by W. Magro et alW^ showed an abrupt transition characterized by a region 
C < for r s =1.86. Later work by B. MilitzeJ^ using more accurate 
variational nodes did not show such a region for r s =1.86. Whether PIMC 
simulations with variational nodes predict a gradual molecular-metallic tran- 
sition for all temperatures, or if the region of an abrupt transition has shifted 
to lower, not yet accessible, temperatures has not been determined. How- 
ever, recent density functional molecular dynamics simulations by various 
authors predict a first order transition below 5 000 

There was a lot of recent interest focused on the hydrogen EOS because 
of the unexpectedly high compressibility inferred from the laser-driven shock 
wave experiment by Da Silva et a/P^ and Collins et aZp^l using the Nova 
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laser at Lawrence Livermore National Laboratory. The measurements indi- 
cated that hydrogen could be compressed to about 6-fold the initial density 
rather than 4-fold as indicated by the Sesame modelP^^ Even though the 
temperatures and pressures reached in shock experiments are higher than 
those in the giant planets these measurements are crucial here since they 
represent the only way to distinguish between different EOS models at high 
temperature. 

The Nova results challenged the existing understanding of high P-T 
hydrogen and triggered many new experimental and theoretical efforts. Dif- 
ferent chemical models gave rise to very different predictionJ^ESEZI ranging 
from 4-fold compression as suggested irP512£l to 6-fold as predicted by the 
Ross modelP^ While the accuracy of chemical models did not allow any 
conclusive predictions, first principle simulations from PIMCP^H as we \\ 
as from density functional molecular dynamics^EEMD consistently predict a 
lower compressibility of about 4.3. 

Since then there have been many attempts to resolve this discrepancy 
but the most important contributions came from new experiments by Knud- 
son et a?! 41 l 42 l 4 ^ at Sandia National Laboratory. Instead of a laser drive, they 
used magnetically driven shock waves in combination with bigger samples. 
They found a significantly lower compressibility quite similar to predictions 
from first principles methods. The new results are also supported by a 
third set of experiments by Russian investigators using spherically converg- 
ing shock waves P=ElD 



4. HYDROGEN-HELIUM MIXTURES 

Due to importance for astrophysical applications, EOS models for hy- 
drogen-helium mixtures have been studied for quite some time. The most 
widely used EOS was derived by Saumon, Chabrier and van Horn (SCH)P^ 
Like previous chemical models, it is based on a hydrogen and a helium EOS 
that combined using an ideal mixing rule. 

Following the discussion of the hydrogen EOS above, we are now ana- 
lyzing the second ingredient: the EOS of helium. Given its low abundance 
in giant planets, helium is not expected to be present in very high concentra- 
tion. However, StevensorP^l has proposed that the hydrogen-helium mixture 
could phase separate under certain high pressure conditions. As a result, 
helium droplets would form and fall as rain towards planet core. The associ- 
ated release of latent heat would delay the cooling process of the planet. This 
was suggested as one possible explanation of why standard models predict 
Saturn to cool at a much faster rate than is observed.^ In order to detect 
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hydrogen-helium phase separation one needs an EOS of pure helium that we 
will now discuss. 




T(1000 K) T(1000 K) 

Fig. 5. Equation of state comparison for pure helium showing the pressure 
(left) and the internal energy (right) for different density as derived from 
path integral Monte Carlo simulations and the chemical model by Saumon, 
Chabrier, and van HornP 



Figure |S] shows an EOS comparison of our PIMC calculations with the 
SCH model. One finds that the pressure is underestimated by SCH and that 
the deviations increase with density. The internal energy at high tempera- 
ture is also underestimated. The analysis suggests a more careful treatment 
of the helium ionization states is needed to improve the chemical model. 

To conclude the comparison with chemical models, we test one more as- 
sumption that is generally made: the ideal mixing hypothesis. We performed 
a large number of PIMC calculations of fully interacting hydrogen-helium 
mixtures at a fixed density of r s = 1.86 for various temperatures and mixing 
ratios. The resulting correction to ideal mixing is given by, 

A/ mix = f{x) - (1 - x) f{x = 0)-xf(x = l). (7) 

Figure HO shows that the corrections to pressure, AP m ; x , increase with 
decreasing temperature reaching 10% for 15 625 K. In the considered temper- 
ature interval, the corrections to the internal energy are largest for 31 250 K 
which can be explained by the different ionization states of the two fluids 
that lead to larger error if ideal mixing is assumed. 

Finally, we discuss pair correlation functions, <?(r), for hydrogen-helium 
mixtures, 
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Helium fraction x Helium fraction x 

Fig. 6. Excess mixing pressure AP m i x and internal energy per electron AE m { x 
are shown for a hydrogen-helium mixture at r s = 1.86. Path integral Monte 
Carlo results are shown for several temperatures. The mixing is performed 
at constant density. A£ m j x is largest for T=31250K because at this tem- 
perature the two end members are characterized by very different degrees of 
ionization. 



g(r) functions are a standard tool to characterize the short-range correlation 
of particles in liquids. In Fig. [JJ we compare the correlation functions for 
pure hydrogen, a x=50% mixture, and pure helium at fixed temperature 
(15 625 K) and density (r s = 1.86). Changes in the pair correlations are 
now discussed as a function of helium concentration. The electron-proton 
g(r) shows a strong peak at the origin due to the Coulomb attraction. The 
peak is not affected when 50% helium is added to a pure hydrogen sample. 
Similarly, the correlation between electrons and helium nuclei is not altered 
by the presence of protons. 

The electron-electron correlation depends strongly on spin. For pairs 
with parallel spin, exchange effects lead to a strong repulsion. The resulting 
g{r) function does not depend much on whether helium is present or not. The 
correlation function between pairs of antiparallel spins, on the other hand, 
is strongly affected by the presence of helium nuclei. In the helium atom, 
two electrons with opposite spin are attracted to the core, which indirectly 
leads to the observed increase in the electron-electron g[r). 

It is interesting to note that proton-proton g{r) changes with the helium 
concentration but the correlation of helium nuclei does not. With increasing 
presence of helium, a peak in the proton-proton g{r) appears at r = 1.4 an 
which indicates the formation of H2 molecules. Adding helium nuclei leads 
to the localization of a fraction of the electrons. The available space in 
combination with the reduced electronic exchange effects then leads to the 
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Fig. 7. Comparison of different pair correlation functions, g(r), from three 
PIMC simulations at r s = 1.86 and 15 625 K: pure hydrogen (red dotted 
lines), a x = 50% mixture (black dashed lines) and pure helium (blue solid 
lines). The six graphs show g(r) for different pairs of protons (p), electrons 
(e) and helium nuclei (a). For electrons, pairs with antiparallel (upper right) 
and with parallel (lower right) spins are distinguished. 

formation of hydrogen molecules which are not present in pure hydrogen at 
the same temperature and density. 



5. CONCLUSIONS 

The discussion in this article centered on the phase diagram of hot dense 
hydrogen. The first PIMC results for hydrogen-helium mixtures were pre- 
sented here, and the accuracy of existing helium EOS models as well as the 
commonly used ideal mixing rule were analyzed. This analysis revealed sub- 
stantial inaccuracies in existing EOS models. These uncertainties prevent us 
from making reliable models for the interior of giant planets and from draw- 
ing conclusions about their formation mechanism, either nucleated capture 
of nebular hydrogerP or gravitational instability driven formation in proto- 
stellar disks.™ 

An improved EOS will also yield a better characterization of the 120 
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extrasolar giant planets that have been discovered so far using radio velocity 
measurements.^ In particular the mass-radius relationship obtained from 
transit extrasolar giant planets^ will allow one to draw conclusions about 
the heavy element abundance in those planets. The metallicity of a planet 
relative to its central star constrains the role of accretion of planetesimals in 
the planet's formation. 
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